% clear workspace
clear

% test data and cov. matrix
mu1 = [5 6 5];
mu2 = [0 1 1];
co = [1 1 1];

% a priori probabilities
prior = [1 1];

% number of data in each set - should match the
% probability distribution above.
anz1=100;
anz2=100;

% generate two datasets
data1 = mvnrnd(mu1, co, anz1);
data2 = mvnrnd(mu2, co, anz2);

% recompute the mean values
mu1 = mean(data1);
mu2 = mean(data2);

% create sample data and group vector
sample = vertcat(data1, data2);
groups = vertcat(repmat(1, anz1, 1), repmat(2, anz2, 1));

% use classify to get the discriminate function using the identity matrix
[C,err,P,logp,coeff] = classify(sample, sample, groups, 'diagquadratic', prior);
% get the coefficients of the discriminate function
K = coeff(1,2).const;
L = coeff(1,2).linear;
Q = coeff(1,2).quadratic;
% build a function for ezmesh
f = @(x,y,z) K + [x y z]*L + sum([x y z] .* ([x y z]*Q), 2);

% visualize
figure
% workaround to enable "hold"
plot3(0,0,0);
hold on

% draw first data cloud
for i=1:anz1
 plot3(data1(i,1), data1(i,2), data1(i,3), 'ro');
end

% draw second data cloud
for i=1:anz2
 plot3(data2(i,1), data2(i,2), data2(i,3), 'bo');
end

% draw the means and the connecting line
plot3(mu1(1,1), mu1(1,2), mu1(1,3), 'gx');
plot3(mu2(1,1), mu2(1,2), mu2(1,3), 'gx');

% concat the two means for line drawing
muline = vertcat(mu1, mu2);
plot3(muline(:,1), muline(:,2), muline(:,3));

% draw the discriminant function
xv = linspace(-5,10,50); % vectors to cover the range of each column
yv = linspace(-5,10,50);
zv = linspace(-5,10,50);
[xx,yy,zz] = meshgrid(xv,yv,zv);
f = @(x,y,z) K + [x y z]*L + sum([x y z] .* ([x y z]*Q), 2);
v = f(xx(:),yy(:),zz(:));
v = reshape(v,size(xx));
isosurface(xx,yy,zz,v,0);

% show grid
grid on
